function [calibrationerror, out] = calibration(searchparameters,fixedparameters,datatargets,gammaf, boundary,penalty)
%This file is the main calibration routine

%Fixed Paramaters =[alpha0;alpha1;alpha2;B;Yss;pi1;gamma;Yss,fringewealth];
B=fixedparameters(1);
Yss=fixedparameters(2);
gamma=fixedparameters(3);
omegaFB = fixedparameters(4);

%Parameters to be calibratged
rho=searchparameters(1);
delta=searchparameters(2);
relativesize=searchparameters(3);
pi1=searchparameters(4);

%Summarize probabilities
pi2=1-pi1;
prob=[pi1,pi2];

%Initial endowment and vector of aggregate endowments
W=Yss;
Yvec=[delta*Yss,Yss];

%vector of bond payoffs
bondpayoff=[rho*B,B]; 

%Back out the income distribution
pAD=prob.*(W./Yvec).^(gamma);%Prices in the first-best Competitive Equilibrium
alpha=omegaFB - bondpayoff./Yvec;%date 1 income shares
alpha0=omegaFB + (pAD*bondpayoff')./W; %date 0 income share
w1=alpha0*W;
w2=W-w1;
y1vec=alpha.*Yvec;
y2vec=Yvec-y1vec;
yf=(Yvec/W).^(gamma/gammaf);

%fsolve optimization parameters
options=optimset('MaxFunEvals',2000,'MaxIter',2000,'TolFun',10^-6,'Display','off');

%%% SOLVE FOR EQUILIBRIUM %%%

%Since we have a good guess for the competitive benchmark, we start at mu=0
%and solve on a grid up to the target relativesize mu/mf.
gridsize=101;
mugrid=linspace(0,relativesize,gridsize);
sol=zeros(gridsize,4);

for ii=1:gridsize
    
    mumf=mugrid(ii);
    
    if ii==1
        guess_complete=[0.1,0.1,0,0];
        
    else
        guess_complete=sol(ii-1,:);
    end
    
    fun_complete=@(controls)eqfind_complete(controls,pi1,w1,w2,y1vec,y2vec,yf,gamma,gammaf,mumf,boundary,penalty);
    [sol(ii,:),~,~]=fsolve(fun_complete,guess_complete,options);
    

    
end

%Compute equilibrium moments
[~,~,bondprice,bidask_bond,~,bidask_cds,~,~,CDS_basis] = output_complete(sol(gridsize,:),pi1,w1,w2,y1vec,y2vec,yf,gamma,gammaf,relativesize,bondpayoff);

bondyield= 100 * (B/bondprice-1);

%Turn bid-ask spreads to basis points
bidask_bond=100 * bidask_bond;
bidask_cds= 100 * bidask_cds;
CDS_basis = 100 * CDS_basis;

%Compute objective 
modelmoments=[bondyield;bidask_bond;bidask_cds;CDS_basis];
diff=modelmoments-datatargets;

calibrationerror= sqrt(diff.^2);

out.alpha = alpha; out.alpha0 = alpha0;

if min(min(alpha),alpha0) < 0 || max(max(alpha),alpha0) > 1
   
calibrationerror = 1e3 * calibrationerror;
    
end

end